Introduction

  • Antimicrobial resistance is becoming a worldwide problem

  • Many pathogens are becoming resistant for many different types of antibiotics

  • The institute of advanced studies (IAS) is doing research to the spreading of gonorrhoea in Amsterdam

  • For gonorrhoeae there might be a new antibiotic available within years

  • It is thus important to find a suitable management plan to prevent new resistance

Mathematical model

  • Aim of our research was to build a network model to model the disease spreading
  • Our network models are based on an existing mathematical model by Fingerhurth et al. (2006) [1]
    • Simple network model
    • Stochastic network

Why use a network model?

  • Sexual network portrays the sexual inter-relationships within a group of people and it is possible to study the spread of STIs in the context of social interactions [2]
  • ODEs, in contrast, fail to capture the complexity of connections between cases
  • By using a network a heterogeneous population can be introduced, which adds more detail in the properties of individuals. This can for example be useful in identifying high-risk groups.
  • Adding this heterogeneity can have influence on the disease dynamics [3]
  • Network properties can for example derived from data, in Amsterdam this is for example possible with GGD data

Mathematical model of antibiotic resistance N. gonorrhoeae spread

In [235]:
Image(filename= "math_model.png")
Out[235]:
In [242]:
plt.figure()
# plt.plot(t, sol_msm[1,:, 0], label='$S_{low}$')
# plt.plot(t, sol_msm[1,:, 1], label='$S_{high}$')
plt.plot(t, sol_hmw[1,:, 2], label='$I_{res, low}$')
plt.plot(t, sol_hmw[1,:, 3], label='$I_{res, high}$')
plt.plot(t, sol_hmw[1,:, 4], label='$I_{sen, low}$')
plt.plot(t, sol_hmw[1,:, 5], label='$I_{sen, high}$')

plt.legend(loc='best')
plt.xlabel('t (years)')
plt.title("Epidemic curve HMW group")
plt.ylabel("Fraction of population")
# plt.ylim(0, 1)
plt.xlim(0, 30)
plt.grid()
plt.show()
In [243]:
plt.figure()
ratio_sol_msm = (sol_msm[:,:, 2] + sol_msm[:,:, 3])/(sol_msm[:,:, 2] + sol_msm[:,:, 3] + sol_msm[:,:, 4] + sol_msm[:,:, 5])
plt.plot(t, ratio_sol_msm[1,:], label = "MSM", c = "g",linewidth=1)
plt.fill_between(t,ratio_sol_msm[0,:], ratio_sol_msm[2,:], alpha=0.5, facecolor = "g")
ratio_sol_hmw = (sol_hmw[:,:, 2] + sol_hmw[:,:, 3])/(sol_hmw[:,:, 2] + sol_hmw[:,:, 3] + sol_hmw[:,:, 4] + sol_hmw[:,:, 5])
plt.plot(t, ratio_sol_hmw[1,:], label = "HMW", c = "b",linewidth=1)
plt.fill_between(t,ratio_sol_hmw[0,:], ratio_sol_hmw[2,:], alpha=0.5, facecolor = "b")
plt.legend(loc='best')

plt.xlabel('t (years)')
plt.ylabel("Fraction of infected individuals")
plt.grid()
plt.ylim(0,1)
plt.show()

Simple Contact network

  • We used a simplified version of the rules in the mathematical model
  • The network is a random erdos renyi network with 500 nodes and an average degree of 10

The parameters

  • Spontaneous healing rate $\nu$
  • Treatment rate $\tau$
  • Hosts develop resistance during treatment with probability $\mu$
  • Transmision probability $\beta$

Build a network

In [264]:
# # create power law graph
# exp = 2.4
# seed = 123456789
# num = 1000
# graph = power_law_graph(exp, num, seed)

# or, create random network
degree = 10
num = 500
p = degree/num
seed = 1234
graph = random_network(p, num, seed)

# compute average degree
degrees = dict(graph.degree())
sum_of_edges = sum(degrees.values())
average_deg = sum_of_edges / len(degrees)
print(average_deg)

# compute r0 (TODO)
9.752
In [265]:
plt.figure()
pos = nx.spring_layout(graph)
nx.draw(graph, pos, node_color='b', node_size=10, with_labels=False)
plt.show()

Implementation of disease spreading

In [268]:
# set person class to each node
for i in range(len(graph)):
    graph.node[i]["Data"] = Individual(i)

# create set of infecteds, initialize with 10 infecteds
init = np.random.choice(len(graph), 10, replace=False)

# time steps
t = 0
t_tot = 20
time_step = 365
dt = 1 / (time_step)

# constants of model
D = np.array([0.19,0.55]) * time_step
mu = 1e-3
PHI = [0.65, 0.50]
BETA = [0.59, 0.87]

# keep track of infecteds and steps
num_steps = int(t_tot / dt)
num_infected = np.zeros((2, num_steps))
num_res = np.zeros((2, num_steps))


p_inf = 0.8
# iterate over sextypes (MSM and HMW)
for sextype in range(2):
    phi = PHI[sextype]
    beta = BETA[sextype]
    nu = (1 - phi)/D[sextype]
    tau = phi/D[sextype]
    beta = 1-(1-1/(average_deg*(1-p_inf)))**(nu)
    num_res[sextype,:], num_infected[sextype,:], graph = network_model(beta, tau, nu, mu, init, num_steps, copy.deepcopy(graph), True)
100%|██████████| 7300/7300 [00:44<00:00, 162.88it/s]
100%|██████████| 7300/7300 [00:42<00:00, 169.92it/s]

Plotting

In [280]:
# plot number of infecteds
t = np.linspace(0,t_tot, num_steps)    
plt.figure()
print("Mean number of infecteds:", np.mean(num_infected))
plt.plot(t, num_infected[0,:], label='all infecteds MSM')
plt.plot(t, num_res[0,:], label='resistant infecteds MSM')
plt.plot(t, num_infected[1,:], label='all infecteds HMW')
plt.xlabel("Time (years)")
plt.ylabel("number of people")
plt.plot(t, num_res[1,:], label='resistant infecteds HMW')
plt.legend()
plt.show()
Mean number of infecteds: 288.4821917808219
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-280-cb51f78d178f> in <module>()
      3 plt.figure()
      4 print("Mean number of infecteds:", np.mean(num_infected))
----> 5 plt.plot(t, num_infected[0,:], label='all infecteds MSM')
      6 plt.plot(t, num_res[0,:], label='resistant infecteds MSM')
      7 plt.plot(t, num_infected[1,:], label='all infecteds HMW')

IndexError: too many indices for array
In [281]:
plt.figure()
plt.plot(t,num_res[0,:]/num_infected[0,:], label='fraction of resistant strain MSM')
plt.plot(t,num_res[1,:]/num_infected[1,:], label='fraction of resistant strain HMW')
# plt.plot(num_res)
plt.xlabel("Time (years)")
plt.ylabel("Fraction of infected individuals")
plt.legend()
plt.show()
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-281-9a4093e5464e> in <module>()
      1 plt.figure()
----> 2 plt.plot(t,num_res[0,:]/num_infected[0,:], label='fraction of resistant strain MSM')
      3 plt.plot(t,num_res[1,:]/num_infected[1,:], label='fraction of resistant strain HMW')
      4 # plt.plot(num_res)
      5 plt.xlabel("Time (years)")

IndexError: too many indices for array

Animation

In [279]:
# set person class to each node
for i in range(len(graph)):
    graph.node[i]["Data"] = Individual(i)

# create set of infecteds, initialize with 10 infecteds
init = np.random.choice(len(graph), 10, replace=False)

# time
t = 0
t_tot = 5
time_step = 365
dt = 1 / (time_step)

# constants of model
D = np.array([0.19,0.55]) * time_step
mu = 1e-3
PHI = [0.65, 0.50]
BETA = [0.59, 0.87]

# choose MSM or HMW to model
sextype = 0
phi = PHI[sextype]
beta = BETA[sextype]
nu = (1 - phi)/D[sextype]
tau = phi/D[sextype]
beta = 1-(1-1/(average_deg*(1-p_inf)))**(nu)

# initial state
_, _, graph = network_model(beta, tau, nu, mu, init, 0, graph, doInit=True, disable_progress=True)

# get values of disease status
values = []
for j in range(len(graph.nodes())):
    if graph.node[j]["Data"].disease_status == 0:
        values.append("Black")
    elif graph.node[j]["Data"].disease_status == 1:
        values.append("Blue")
    elif graph.node[j]["Data"].disease_status == 2:
        values.append("Red")
        
# create figure and set axes
fig = plt.figure(figsize = (10,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
im1, = ax2.plot([], [], c = 'b', label='all infecteds MSM')
im2, = ax2.plot([], [], c = 'r', label='resistant infecteds MSM')
ax2.legend()
ax2.set_xlim(0,t_tot)
ax2.set_ylim(0,num)

# initialize
pos = nx.spring_layout(graph) 
plt.close()
ani_step = 365//5
frames = (t_tot*365)//ani_step
num_infected = np.zeros(ani_step*frames)
num_res = np.zeros(ani_step*frames)

# update function for animation
def update(i):
    if i<1:
        temp_res, temp_infected, _ = network_model(beta, tau, nu, mu, init, 1, graph, doInit=False, disable_progress=True)
    else:
        temp_res, temp_infected, _ = network_model(beta, tau, nu, mu, init, ani_step, graph, doInit=False, disable_progress=True)
    num_infected[i*ani_step:(i+1)*ani_step] = temp_infected
    num_res[i*ani_step:(i+1)*ani_step] = temp_res
    t = np.linspace(0,i*ani_step/365, i*ani_step) 
    im1.set_data(t, num_infected[:i*ani_step])
    im2.set_data(t, num_res[:i*ani_step])
    ax1.clear()
    
    # solution at step i
    values = []
    for j in range(len(graph.nodes())):
        if graph.node[j]["Data"].disease_status == 0:
            values.append("Black")
        elif graph.node[j]["Data"].disease_status == 1:
            values.append("Blue")
        elif graph.node[j]["Data"].disease_status == 2:
            values.append("Red")
    nx.draw(graph, pos=pos, node_color=values, node_size=10, ax=ax1, clim = 2)
    
    # Scale plot ax
    ax1.set_title("Frame %d:    "%(i+1), fontweight="bold")
    ax1.set_xticks([])
    ax1.set_yticks([])
    
# animate!
animation.FuncAnimation(fig, update, frames=frames, interval=100, repeat=True)
25
Out[279]:


Once Loop Reflect

Stochastic network model

The stochastic network model is based on an existing model by Kretzschmar, M. et al. (1996) [4]

The system consists of Persons defined by:

  • Gender (male / female)
  • Age (15 - 64)
  • Sexual activity (high / low)
  • Disease status (uninfected or (non-)resistant (a)symptomatic)

During each timestep, Partnerships are formed between two Persons. A Partnership is defined by:

  • Type (casual or steady)
  • The two Persons involved

The sexual contact network is therefore dynamic in time

Each timestep of the model consists of seven parts:

  1. Partnership formation
  2. Disease transmission
  3. Separation of partnerships
  4. Replacement
  5. Recovery
  6. Treatment
  7. Advance time

Initialize model

  1. Initialize 10,000 Persons of random age. 5% of Persons with age < 35 have high sexual activity
  2. Perform 2,000 iterations to obtain the steady state
  3. Infect 10% of the core group (high sexual activity) with non-resistant asymptomatic disease

The resistant disease enters the system through treatment; treatment will result in contracting the resistant disease with probability 0.01%

Time-evolution towards the steady state

In [231]:
fig

Time-evolution of the population

In [ ]:
fig

Time-evolution of infecteds over time

In [ ]:
fig

Fraction of resistant infecteds over all infecteds over time

In [ ]:
fig

Infectiousness of the disease ($R_0$)

Although the disease might seem quite infectious, simulations have shown $R_0=0.90 \pm 2.93$ (N=1,000, 95% CI). The spread of the disease is very dependent on where the disease is started. If a person with high sexual activity is infected, it will quickly spread throughout the entire core group whereas a person with low sexual activity may only affect a select few different partners

Animation

In [211]:
anim
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-211-027b5d0310db> in <module>()
----> 1 anim

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/displayhook.py in __call__(self, result)
    255             self.start_displayhook()
    256             self.write_output_prompt()
--> 257             format_dict, md_dict = self.compute_format_data(result)
    258             self.update_user_ns(result)
    259             self.fill_exec_result(result)

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/displayhook.py in compute_format_data(self, result)
    149 
    150         """
--> 151         return self.shell.display_formatter.format(result)
    152 
    153     # This can be set to True by the write_output_prompt method in a subclass

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/formatters.py in format(self, obj, include, exclude)
    178             md = None
    179             try:
--> 180                 data = formatter(obj)
    181             except:
    182                 # FIXME: log the exception

<decorator-gen-9> in __call__(self, obj)

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/formatters.py in catch_format_error(method, self, *args, **kwargs)
    222     """show traceback on failed format call"""
    223     try:
--> 224         r = method(self, *args, **kwargs)
    225     except NotImplementedError:
    226         # don't warn on NotImplementedErrors

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/formatters.py in __call__(self, obj)
    343             method = get_real_method(obj, self.print_method)
    344             if method is not None:
--> 345                 return method()
    346             return None
    347         else:

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in _repr_html_(self)
   1484             return self.to_html5_video()
   1485         elif fmt == 'jshtml':
-> 1486             return self.to_jshtml()
   1487 
   1488 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in to_jshtml(self, fps, embed_frames, default_mode)
   1467                 self.save(f.name, writer=HTMLWriter(fps=fps,
   1468                                                     embed_frames=embed_frames,
-> 1469                                                     default_mode=default_mode))
   1470             # Re-open and get content
   1471             with open(f.name) as fobj:

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs)
   1257                     for anim, d in zip(all_anim, data):
   1258                         # TODO: See if turning off blit is really necessary
-> 1259                         anim._draw_next_frame(d, blit=False)
   1260                     writer.grab_frame(**savefig_kwargs)
   1261 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in _draw_next_frame(self, framedata, blit)
   1294         # post- draw, as well as the drawing of the frame itself.
   1295         self._pre_draw(framedata, blit)
-> 1296         self._draw_frame(framedata)
   1297         self._post_draw(framedata, blit)
   1298 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in _draw_frame(self, framedata)
   1812         # Call the func with framedata and args. If blitting is desired,
   1813         # func needs to return a sequence of any artists that were modified.
-> 1814         self._drawn_artists = self._func(framedata, *self._args)
   1815         if self._blit:
   1816             if self._drawn_artists is None:

<ipython-input-210-e04a745e332f> in plot_graph(i)
     29     if i > 0:
     30         for j in range(365):
---> 31             s.time_step()
     32 
     33     persons = set()

~/Documents/GitHub/CSS2018/Code/stochastic_network_model.py in time_step(self)
    645                 if np.random.random() < self.recovery_rate_asymptomatic_men:
    646                     p.cure(self)
--> 647             elif p.gender == 1 and p.disease_status in [1,3]:
    648                 if np.random.random() < self.recovery_rate_symptomatic_women:
    649                     p.cure(self)

KeyboardInterrupt: 

References

[1] Fingerhuth, S. M., Bonhoeffer, S., Low, N., & Althaus, C. L. (2016). Antibiotic-resistant Neisseria gonorrhoeae spread faster with more treatment, not more sexual partners. PLoS pathogens, 12(5), e1005611.

[2] De, P., Singh, A. E., Wong, T., Yacoub, W., & Jolly, A. M. (2004). Sexual network analysis of a gonorrhoea outbreak. Sexually transmitted infections, 80(4), 280-285.

[3] Bansal, S., Grenfell, B. T., & Meyers, L. A. (2007). When individual behaviour matters: homogeneous and network models in epidemiology. Journal of the Royal Society Interface, 4(16), 879-891.

[4] Kretzschmar, M., van Duynhoven, Y. T., & Severijnen, A. J. (1996). Modeling prevention strategies for gonorrhea and chlamydia using stochastic network simulations. American Journal of Epidemiology, 144(3), 306-317.